
\u 

& GENERALIZED TENSOR FORMULATION 
OF ATMOSPHERE AND SEAS DYNAMICSj 


Lee M. Avis, Richard E. Turner, and Thomas H. Rees 


Langley Research Center 
Hampton, Va. 23665 

i " 

NATIONAL AERONAUTICS AND SPACE ADMIN^pkATiTN 


WASHINGTON, D. C. 



• |SEP££#»i**r75 


TECH LIBRARY KAFB, NM 



TECH LIBRARY KAFB, NM 




OObBM 


i. 


4. 


7. 


9. 


Report No. 

NASA TR R-441 

Title and Subtitle 


2 . 


Government Accession No. 


3. Recipient's Catalog No. 


5. Report Date 


A GENERALIZED TENSOR FORMULATION OF 
ATMOSPHERE AND SEAS DYNAMICS 

Author(s) 

Lee M. Avis, Richard E. Turner, and Thomas H. Rees 

Performing Organization Name and Address 


September 1975 

6. Performing Organization Code 


8. Performing Organization Report No. 

L-9928 


10. Work Unit No. 

176-24-31-02 


NASA Langley Research Center 
Hampton, Va. 23665 


12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D.C. 20546 

15. Supplementary Notes 


16. Abstract 


11. Contract or Grant No. 


13. Type of Report and Period Covered 

Technical Report 

14. Sponsoring Agency Code 


A new generalized mathematical model has been developed for simulation of the 
dynamics and transport of both the atmosphere and seas. A nearly horizontal bottom 
coordinate surface conforms to the land -to -air interface and the sea-floor-to-water inter- 
face, which simplifies computations. General vertical motion of the other quasi-horizontal 
coordinate surfaces is allowed; thus, external gravity waves can be represented by the top 
coordinate surface, and meteorological fronts and inversion layers in the atmosphere and 
refractive layers in the seas can be represented with enhanced resolution by the internal 
quasi -horizontal coordinate surfaces. A tensor reformulation of the standard subgrid 
mixing theory departs significantly from the standard theory in allowing, as a solution under 
adiabatic conditions, rigid -body rotation of the atmosphere. 


17. Key Words (Suggested by Author(s) ) 
Fluid mechanics 

Time -dependent metrics 
Geophysical fluids 


18. Distribution Statement 

Unclassified - Unlimited 


Numerical simulation 


Subject Category 45 


19. Security Cl assif. (of this report) j 

20. Security Classif. (of this page) 

21. No. of Pages 

22. Price* 

Unclassified 

Unclassified 

64 

$4.25 


For sale by the National Technical Information Service, Springfield, Virginia 22161 



A GENERALIZED TENSOR FORMULATION OF 
ATMOSPHERE AND SEAS DYNAMICS 

Lee M. Avis, Richard E. Turner, 
and Thomas H. Rees 
Langley Research Center 

SUMMARY 

A new generalized mathematical model has been developed which is applicable to 
simulation of the dynamics and transport of both the atmosphere and seas. Considerable 
conceptual and computational simplification is achieved by requiring the nearly horizontal 
bottom coordinate surface to follow the land-to-air interface and the sea-floor -to -water 
interface. The model allows general vertical motion of the other quasi -horizontal coordi- 
nate surfaces; thus, external gravity waves can be simulated by letting the top coordinate 
surface ride with a top free surface. Also, the freedom of motion of the internal quasi- 
horizontal surfaces provides a potential for tracking meteorological fronts and inversion 
layers in the atmosphere and refractive layers in the seas with resolutions superior to 
those realizable by simulation models having fixed coordinate lines, for equal "horizontal" 
grid spacing and an equal number of "horizontal" surfaces. 

The standard subgrid mixing theory used by meteorologists is reformulated in terms 
of a strain-rate tensor for both three-dimensional and hydrostatic quasi-two-dimensional 
dynamics. The tensor reformulation departs significantly from the standard theory in 
allowing, as a solution under adiabatic conditions, rigid-body rotation of the atmosphere. 

INTRODUCTION 

There are many advantages in developing computer models of the Earth’s atmos- 
phere or seas in coordinate systems allowing generally specified vertical motion of the 
coordinate gridpoints. In the atmospheric troposphere, coordinate gridpoints can be 
moved vertically to follow the motions of sharp temperature gradients occurring in 
frontal systems to obtain improved resolution of the frontal system. In regions where 
low-pressure cells penetrate into the stratosphere, the coordinate gridpoints can be 
moved vertically to follow the zone separating the turbulent troposphere from the 
thermally inverted stratosphere and thereby reduce false numerical diffusion between the 
two regions. The seas also have a time -dependent vertical structure that can be well 
represented by coordinate systems allowing vertical motion of the coordinate gridpoints. 



In coastal regions where fresh -water rivers empty into the seas, strong vertical 
density stratifications occur and cause time-varying changes in the flow fields. Spurious 
numerical diffusion through such strong vertical density gradients can be reduced by 
moving gridpoints vertically to follow the density gradients. In the case of coastal-zone 
tidal waves, tidal currents can be simulated in conjunction with moving the upper level 
gridpoints vertically to follow the water’s free surface. Vertical density stratifications 
occur along the common boundaries of adjacent oceans (the Arctic and Atlantic Oceans, 
for example), and the open seas are generally stratified because they are top heated. 

The National Aeronautics and Space Administration is presently involved in 
developing computer models of the Earth’s atmosphere and seas to aid in determining 
satellite sensor requirements and to aid in the interpretation of remotely sensed data in 
connection with environmental pollution problems. Since pollution occurs over a wide 
range of distance scales in the atmosphere and seas, a mathematical model of geophysical 
circulation allowing vertical motion of coordinate gridpoints and having validity over a 
wide range of scales is needed for the development of circulation models. 

The literature shows that a number of successful attempts have been made to model 
the atmosphere and seas. References 1 to 4 are descriptions of global atmospheric 
models; reference 2 also describes a global seas model. A survey of small-scale models 
is given in reference 5. 

Generally, atmospheric-circulation models have been in existence longer than ocean- 
circulation models and, consequently, have reached a higher level of sophistication. The 
atmospheric models given in references 2 and 4 use the ”sigma coordinate” method where 
vertical motion of the gridpoints is used to keep constant mass ratios in the model vertical 
layers. Such a model cannot be converted easily to track surfaces of constant tempera- 
ture, for example, and even though the TT sigma coordinate” method allows vertical motion 
of the gridpoints, it cannot be considered a flexible technique. 

The coordinate system and equations of motion of seas models (e.g., see ref. 2) 
generally differ considerably from the coordinate system and equations of motion of 
atmospheric-circulation models (e.g., see refs. 2 and 4). This fact probably increases 
the difficulty of communication between atmospheric and ocean researchers. 

In this report an attempt is made to develop the equations of motion in a coordinate 
system that allows generally specified vertical motion of the gridpoints. The resulting 
equations of motion are cast in terms of dependent variables which represent either 
atmospheric or ocean circulation. Computer models developed from such a formulation 
need differ only in respect to the physics to be included in each model. Consequently, 
the terminology gap between atmospheric and oceanic researchers is bridged. 
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Present-day understanding of mixing processes is poor, and it is anticipated that 
turbulent mixing theories will become considerably more sophisticated, especially for 
pollution transport. Tensor calculus (having been developed from principles of invariance) 
is a good candidate for the description of better mixing theories and is employed in this 
report. 


SYMBOLS 

A determinant of the quasi -horizontal metric tensor 

a determinant of the metric tensor 

(a c ). coriolis and centripetal acceleration along dx 1 

C v heat capacity per mole at constant volume 

Drp,Dg elements in the subgrid diffusion model, defined by equations (D9) 

E internal energy per unit mass 

g acceleration due to gravity 

gj gravity component along dx 1 

gjjjg 1 ’* covariant and contravariant elements of the metric tensor, respectively 

h physical height above reference level of parametric (grid) lines 

Kjj horizontal kinematic viscosity 

Von Karman constant 

Ky vertical kinematic eddy viscosity for momentum 

Ljj physical components of the deviatoric strain-rate tensor 



mixed tensor components of the deviatoric strain-rate tensor 
mixed tensor components of the strain-rate tensor 
unit vector normal to dS, positive outward 
hydrostatic pressure 

general quantity which assumes the role of the dependent variables 

heating rate in the thermodynamic equation (E5) 

water-vapor mass fraction (specific humidity) 

universal gas constant 

Earth's radius 

area 

quasi-source term for transport equation ^for a scalar Q, 
s(^pQ) = S'(^pQ)) 

source term for transport equation 

physical vertical velocity of coordinate gridpoints relative to rotating Earth 
contravariant components of gridpoint velocities relative to rotating Earth 
temperature, °K 

general contravariant and covariant vectors, respectively 
time, sec 

physical components of velocity vector relative to rotating Earth 


physical components of turbulent velocity vector relative to rotating Earth 



uj covariant tensor component of velocity vector relative to rotating Earth 

Vj physical components of fluid velocity relative to coordinate frame of reference 

v* contravariant form of Vj 

X 3 depth of atmosphere or seas in transformed coordinates 

1 2 

x 1 coordinate variables for reference coordinate system, where x A and x 

lie along the real Earth surface and x 3 is perpendicular to the real 
Earth surface 

z^ physical coordinate in direction of x 1 

at, /3 indices restricted to the horizontal surface 

y ratio of heat capacity at constant pressure to heat capacity at constant volume 

A coordinate increment 

6 physical height of Earth's terrain from a given reference level 

6 characteristic 6 on a grid interval 

sj mixed tensor components of Kronecker delta 

error fraction in the square of a line element 
£ generalized density defined by p 

1 0 Q 

dy. coordinate volume element, dx dx dx 

£ variable of integration 

p density of air or water 

a > Earth's spin rate 

uT mean molecular weight 
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Operators: 


space -time average over a coordinate volume v and time increment At, 

-L_ f f ( ) dx 1 dx 2 dx 3 dt 
v At J^t *^Af 

denotes turbulent part of ( ) 
square of ( ) 

partial derivative of ( ) with respect to x 3 
covariant derivative of ( ) with respect to x^ 

co-moving time derivative of ( ) 

partial derivative of ( ) with respect to t 

Christoffel symbols of the second kind 

ANALYSIS 

Introductory Remarks 

A simple, generally valid mathematical model for simulating atmospheric and 
oceanic dynamics and pollution transport was developed with a general, time -varying 
vertical grid structure. A review of the literature (e.g., refs. 1 to 4, 6, and 7) indicated 
that no known existing published model with generally prescribed time-varying vertical 
grid structure had been applied to both atmospheric and oceanic circulation. Reference 1 
describes an atmospheric-circulation model wherein the vertical coordinate is the fixed 
physical distance from sea level. References 2, 3, and 4 describe atmospheric-circulation 
models in the "sigma coordinate" wherein the ratio of pressure to surface pressure is 
used as the vertical coordinate. The "sigma coordinate" technique allows vertical motion 
of the gridpoints; unfortunately, the vertical motion cannot be prescribed generally. Ref- 
erence 6 describes a coordinate system wherein general vertical motion of the gridpoints 
is allowed; however, the resulting transport equations contain elements from two coordi- 
nate systems and may present difficulties when modeling subgrid transport. The strong- 
point of the approach given in reference 6 is that it is not limited by the small slope 
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approximation for the horizontal coordinate surfaces. Reference 7 describes the "sigma 
coordinate" approach in detail. None of the references gave formulations intended to 
cover both atmospheric and oceanic circulation. 

Some notable features of the new mathematical model are as follows: 

(1) Combining density and the vertical scale factor to form a generalized density 

so that the atmosphere and seas have the same motion equations 

(2) Distortion of the quasi-horizontal bottom coordinate surface to follow the real 

Earth terrain and provision for general vertical motion of the other quasi- 
horizontal coordinate surfaces 

(3) Projection of gravity onto the quasi -horizontal coordinate surfaces 

(The gravity force is derived from the gravity potential specified along 
the quasi -horizontal parametric lines.) 

(4) The representation of the intuitive subgrid diffusion (see ref. 3) as a true tensor, 

which corrects a deficiency in the subgrid diffusion formulation in allowing 
rigid-body rotation of the atmosphere under adiabatic conditions 

(5) The extension of the atmospheric top to zero pressure as in the "sigma 

coordinate" approach (ref. 7), which has the ratio of pressure to surface 
pressure as the vertical coordinate (see also ref. 6) 

(6) The use of generalized density rather than pressure as a dependent variable 

The development of the mathematical model in tensor notation is accomplished in 
appendixes A to F. A brief outline of the development sequence in the appendixes which 
will serve to give them continuity is given as follows: 

In appendix A, the transport equations for momentum and a general motion integral 
are derived from conservation integrals in which the volume elements are allowed to vary 
with time: that is, the coordinate points are allowed to move in some arbitrarily speci- 
fied manner. As a direct result, the motion equations can simulate any conceivable mix- 
ture of Eulerian and Lagrangian viewpoints. 

In appendix B, differential distances for real Earth topography are shown to be 
quadratically dependent on the slope of the real Earth surface. Thus, for small surface 
slopes the geometry is, to second order, unchanged by mountains and sea bottoms, and, 
to the order of the square of the slope of the real Earth geometry, the equations of motion 
for geophysical circulation depend on the real Earth topography only through the external 
forces, such as gravity. 

In appendix C, one transport equation, valid for all motion integrals, including sub- 
grid transport terms and appropriate source terms, is given. 
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In appendix D, the subgrid transport treatment of Washington and Kasahara (ref. 1) 
is converted into tensor notation for full three-dimensional turbulence transport and is 
then simplified to the usual two-dimensional stratified turbulence transport. The tensor 
elements of the strain-rate tensor differ from the strain-rate matrix of reference 1 with 
the result that the tensor treatment of subgrid momentum transport allows rigid-body 
rotation of the atmosphere under adiabatic conditions. An alternate tensor treatment of 
subgrid momentum transport, also allowing adiabatic rigid-body rotation, is reported in 
reference 8. 

In appendix E, a tensor mathematical model for full three-dimensional dynamics 
with full three-dimensional Reynolds transport in physical components, as opposed to 
tensor components, is presented. 

Finally, in appendix F, the usual hydrostatic assumption is applied to the momentum 
equations. (See ref. 1, for example.) As a direct result of the hydrostatic assumption, 
one can calculate pressure from the known density field. The hydrostatic assumption 
causes the atmosphere and oceans to be similar except in respect to computation of the 
vertical velocities. 

In the nearly incompressible oceans, the vertical velocity v^ may be computed 
from the mass-conservation equation, and because (see ref. 9) work of compression is 
safely ignored in the internal -energy balance, one need not compute the divergence of the 
velocity uh 

Dynamic simulations of the atmosphere generally include work of compression in 
the internal-energy balance. If one chooses pressure rather than density as the dependent 
variable to be simulated, then the internal-energy constraint is automatically satisfied; 
however, then one must be careful to conserve mass. In the present report, density is 
chosen over pressure as the dependent variable and, consequently, mass is automatically 
conserved, whereas the internal -energy constraint is satisfied only if one computes the 
vertical velocity from a conservative form of the internal -energy equation. 

Conservation Equations for the Atmosphere and Seas 

Over a wide range of scales, the atmosphere and seas approximately satisfy the 
hydrostatic assumption, for which condition the atmosphere and seas have a common set 
of conservation equations for mass, momentum, and energy. The atmosphere has water 
vapor as an additional variable, whereas the seas have salt. 

The coordinate system chosen to present the conservation equations is an 
approximately orthogonal coordinate set wherein the coordinate surfaces defined by 
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X 3 = Constant are nearly spherical, with the x 3 = 0 coordinate surface following the 
real Earth terrain. All other "horizontal" parametric surfaces are allowed to be posi- 
tioned quite generally, or even to be in motion, but with the restriction that the slopes of 
all "horizontal" coordinate surfaces are not very much greater than the maximum slope 
of the real Earth terrain. The x* coordinate is the co-latitude angle, the x 3 coordi- 
nate is the longitude angle, and the x 3 coordinate direction is locally perpendicular to 
the quasi-horizontal coordinate surfaces. 

In appendix B, it is shown that the differential geometry of the chosen coordinate 
system is sufficiently represented by the approximate metric tensor 

g n = (r ) 2 

g 22 = ( r sin x 1 ) 

§33 = § 33 ( x1 ’ x2 ’ x3 ’ 0 


where r is the Earth's radius; and, thus, the determinate of the quasi -horizontal metric 
tensor becomes 


A " g ll g 22 


From appendix E, one has a common set of conservation equations for the atmos- 
phere and seas wherein the real density p is combined with the vertical scale factor 
^g 3 2 to form a generalized density £ = p^ggg. 

The resultant mass-conservation equation in terms of physical velocities Vj is, 
from equation (E2), 



4=-yi/Ac 

S/A 8x a 




= 0 


The momentum equation along the x 1 parametric line is, from equations (E19) 
and (E21), 


( 1 ) 
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= /g 33 3P C TJ 9S 3 j.^1 r g 9h 

' ^I 9 - 1 39x1 V^T v^ 9xl 

1 a U (VlVg)\ B / 

^ 9X "\ / 9x3 \ \fih$3 


- Cl 


< U 2’ U 2 , > 


cxcx 

U„ 


+ V, 


cot X* 


( 2 ) 


The momentum equation along the x^ parametric line is, from equations (E22) and (E24), 

1 9 / /— Vn, \ o / Vc 


= _\/HI iH + c y 8S 3 , ^2 _ , g 3h 


^a x 2 ^ 'ax 2 


3x^ 


i _»_Wa C 

(Ui’u 2 '> u. 


< U a ’ U 2')\ 8 4 < U 2’ U 3') 


■75C 

ax" 


g 33 


- c 


+ Vi \ 1 cot x* 


fa 'fa i 


(3) 


where the Greek indices take on the values 1 and 2 and, from reference 1, 

( a c)l ! 

— - - = -2o>U 9 cos x 

fa 


( a c) 


V s 22 


2 1 
— = 2u>U^ cos x 


(4) 
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and, it should be noted, where the centrifugal term proportional to (u>) has been con- 
sidered to be a part of the gravity term, which is the usual practice. 

The internal -energy conservation equation is, from equation (E5), 



The final conservation equation for atmospheric water or ocean salt, or for an arbitrary 
pollutant, is, from equation (E6), 



( 6 ) 


where 



is the source term for water vapor resulting from evaporation con- 


densation. For the seas, the source term would reflect salt dissociation precipitation; 
for pollution, the source term would represent scavenging and chemical reactions. 


Subgrid Correlations 

The subgrid mixing model is the tensor representation of an intuitive subgrid model 
developed by Smagorinsky and described in reference 3. Simply put, subgrid waves, or 
waves unresolvable on the grid spacing, have a sizable effect on the development of circu- 
lation patterns. It is implicit in Smagorinsky' s model that these subgrid waves act like 
diffusion and tend to drive property gradients to zero and, therefore, represent the irre- 
versible nature of circulating fluids. The horizontal subgrid action is similar to wind- 
or current-driven turbulence, whereas the vertical subgrid action is similar to thermally 
driven turbulence. 
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The horizontal subgrid mixing, then, is taken to be a function of the strain-rate 
tensor, which is a measure of the mechanical stirring of the fluid by the resolvable waves. 

The vertical mixing in the atmosphere (being driven by small-scale convective over- 
turning, as, for example, by cumulus convection) is taken to be a function of the thermal 
lapse rate. (See ref. 1.) The vertical mixing the oceans (which is top heated and, there- 
fore, generally stable thermally) is modeled more like mechanical turbulence in refer- 
ence 9 by specifying a very low eddy viscosity in the deep ocean and specifying relatively 
high eddy viscosity in the surface wind-driven currents where mechanical stirring is 
strong. 

The subgrid mixing model begins with the deviatoric strain-rate tensor which is a 
measure of the stirring of the fluid. The deviatoric strain-rate tensor 1 j from appen- 
dix D is, for the quasi-two-dimensional case, 
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Horizontal Subgrid Mixing 

The strain-rate tensor is related to the subgrid horizontal-velocity correlation 
tensor as in reference 3 by 


<W) - -Wap “» 

where represents the physical components of the strain-rate tensor and is 

computed by 


^a/3 ^./3 ^ 

One simple interpretation of equation (8) is that velocity gradients in the fluid (measured 
by cause a potential for momentum to be transferred by means of the velocity 

correlation (u a 'u 3 ') , and K h is simply a constant of proportionality. The horizontal 
subgrid mixing for a scalar Q, as in reference 1, is 


(V 6 *’) 


k h aQ 



( 10 ) 


where gradients in Q cause a potential for Q to be transferred and, again, Kjj is the 
constant of proportionality. 

The kinematic eddy viscosity is related to the strain-rate tensor for mechanical 
stirring by 


K H = '^( K o a ) 2 (‘^ J ! 3 a,) 1/2 <“> 

where Kq is the Von Karman constant and A is the grid spacing. Equation (11) also 

has a simple interpretation. Mechanical stirring ^measured by \j f causes the 

eddy viscosity, whereas ^ is a constant of proportionality which allows one to 

time his subgrid mixing model to best simulate the real world. Equation (11) implies 
that, for simulation cases where A approaches zero, all important mixing processes 
can be resolved on the grid. 
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Vertical Subgrid Mixing 


In the case of vertical transport, the actual mixing process is generally thermally 
driven in the atmosphere. Reference 1 gives a parameterized description of the vertical 
kinematic eddy viscosity Ky, which is dependent on the temperature lapse rate. In the 
oceans, Cox (ref. 9) has two values of Ky. In the deep ocean where mechanical stirring 
is low, a small value of Ky is used; whereas for the upper ocean in the wind-driven 
currents where mechanical stirring is strong, a large value of Ky is used. 

The deviatoric strain-rate tensor is related to the subgrid velocity correlations 
corresponding to vertical mixing of the horizontal momentum, as in reference 1, by 


<V U 3’> - - K V £ c3 


( 12 ) 


where, again, 


L 


a3 




(13) 


and equation (12) has the same interpretation as equation (10). The vertical subgrid 
mixing for a scalar Q, as in reference 1, is 



K V 9Q 

^ 9 ? 


(14) 


Vertical-Velocity Computations 


The prime unknown remaining to be specified is the vertical velocity. In hydro- 
static models such as those given herein, the vertical velocity cannot be computed from 
the vertical-momentum equation. 


p i 

Seas .- In ocean-circulation models the work of compression — ujj ■ is customarily 

neglected in the internal -energy conservation statement given in equation (5) which 
implies that 


8U3 

9X 3 



(15) 
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Clearly, at the bottom (x 3 = o), the vertical velocity Ug vanishes, and so 


V3= ~hh 


u. 


\J & aal 




( 16 ) 


where | is a variable of integration. The vertical velocity Ug is not required to solve 
the ocean-circulation equations (as approximated herein); however, values of Ug might 
be a desired output of the simulation. One has the freedom to assign Vg or Sg for 
input to the transport equations, and for computing £. 

An additional approximation is employed, the Boussinesq approximation, which states 
that density variations are negligible except when multiplied by the acceleration of gravity. 

Atmosphere.- In atmospheric simulations, the work of compression is usually con- 
sidered in the internal-energy equation, and so computation of Ug is required for solu- 
tion of the equation for atmospheric models; however, when Ug has been computed, one 
still has the freedom to assign values of either Vg or Sg and, consequently, move the 
vertical position of the horizontal coordinate surfaces to any desired position. 


Several thoughts which should be considered before the vertical -velocity computa- 
tional sequence is developed are given as follows: 

(1) The pressure can be computed from the density profile by utilizing the hydro- 
static assumption. 


(2) The internal energy can be computed from pressure and density from the equa- 
tion of state. 


(3) Since internal energy can be computed independent of equation (5), and since 
equation (5) must be satisfied for consistency of internal energy, then equation (5) should 
be used to compute Ug. 

(4) The physical velocity component Vg can then be chosen arbitrarily to give any 
desired vertical coordinate spacing; or any desired vertical distribution of £ can be 

f x3 3 

chosen equally well provided \ £ dx J is properly conserved, because this integral is 

the total mass per unit cross section of area along the nearly vertical x 3 coordinate 
line. With these thoughts in mind, one can proceed generally to specify Vg and then 
compute 9£/8t and Ug; or one can specify 8£/8t and then compute Vg and Ug 
equally well. In either case, the boundary conditions on Vg at the bottom (x 3 = 0) and 
the top (x 3 = X 3 ) for a conservative -mass system are 

V 3 (°) = V 3 (x 3 ) = 0 (17) 
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In illustration of these considerations, Vg is chosen, subject to the boundary con- 
ditions given in equation (17). Next, 9£/9t is computed from the mass -conservation 
equation (1). The final task is to compute Ug from equation (5). In order to solve 

equation (5) for Ug, one must know (£E) and ^~(£E). 

The computation of (£E) begins with the hydrostatic assumption where 



(18) 


Then, from the equation of state, 


p _ PRT _ ?RT 

" V*3S“ 


(19) 


and since 


E = 



( 20 ) 


equations (19) and (20) combine to give 



( 21 ) 


Consequently, the hydrostatic approximation allows one to convert equation (5) from (£E) 
description to P dependence, and thereby simplify the solution, because P and 8P/9t 
can be computed from simple density integrals. The right-hand side of equation (5), 


being already in its proper form, can be denoted by 

l/A 

tion (21) into the left-hand side of equation (5) and get 



; then, substitute equa- 



The next step is to get an expression for 9P/9t from equations (1) and (18) which 
results in 
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The first term in equation (22) expands to 



9P 

at 


Then, from equation (Ell), 


Next, substitute equations (23) and (25) into equation (24) and get 

§ +g?v 


_a_ 

at! 


Finally, substitute equation (26) into equation (22), expand — to 

av c 


Vo + pZl 
3 3 ax3 


ax 


set 


and 


= -g? 

ax’* 


V 3+ S 3 = u 3 


and then get 


v 

R 


/M f r* 3 W,/x r -Z° l\ df ^.p!^ i -2-< f ,/A^ P 

v/ A Jy3 Ay® l /cr J 3X^ \/A 3X^ ' ^ 


v/A ^x^ ax' 

i a/.fr N 


\ 



Equation (27) now involves Ug as the only unknown; thus, replace the right-hand side of 
equation (5) for -y— S^\/A £e), solve for 8Ug/3x^, and get 



IP 


au. 


ax 1 ’ 


c v^«C x3 a /rf, 

-R ^ Jx3 8x“^ 




v a\ P 8 /nr-.rr u « 




V 




aai 


+ CQ - 




< EU 3>\ 
\/^33 / 


(28) 


Equation (28) is then integrated to get Ug by use of the boundary condition, Ug = 0 
at x^ = 0. Equation (28) was derived from the conservation equation for internal energy; 
it is a form of the Richardson equation in reference 1. 

A brief recapitulation of the chosen approach to computing the vertical velocity is 
given as follows: First, Vg, the vertical fluid velocity relative to the coordinate system, 

is specified (consistent with the boundary conditions of eq. (17)). Second, 3?/at is com- 
puted from equation (1). Third, U 3 is computed from equation (28). When V 3 is 
chosen arbitrarily, \/ggg will usually be time dependent, the variation of which can now 
be computed from equation (25) by 


^)=i( S 3)=^( U 3-^) 


(29) 


There are many other valid approaches for specifying the vertical structure of an atmos- 
pheric model, but the important point to be made is that one can choose the vertical struc- 
dure in any desired manner. Equation (28) is the general link by which the vertical struc- 
ture is made to evolve with internal energy properly conserved. 

General Solution Considerations 

The solution to a given circulation problem is typically straightforward in concept. 
The conservation equations describe the time evolution of the chosen dependent variables 
by means of dynamic equations driven by the dependent variables and subgrid correlations 
between the dependent variables. In the present model, approximations have been made 
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on the terms containing real world terrain to gain simplicity in the governing equations 
while incurring only insignificant errors. The results indicate that the real world is, 
for practical grid spacings, equivalent to a spherical world where the gravity vector is 
oblique to the spherical surface. The resulting model allows one to specify general ver- 
tical grid structures for the model as a function of time. The model was developed 
specifically to simulate both atmosphere and seas motions. 

CONCLUDING REMARKS 

A generally valid, but simple, mathematical model suitable for simulating atmos- 
pheric and oceanic dynamic motions and pollution transport is presented. The mathe- 
matical model was formulated in general tensor notation, and in the conservation equa- 
tions, the density p and vertical scale factor y/gjjjj 1 have been grouped together to 
form a generalized density variable ? suitable for oceanic or atmospheric simulation. 
For simplifying computations the quasi-horizontal bottom coordinate surface has been 
deformed to follow real Earth terrain, whereas the top coordinate surface moves verti- 
cally with the oceanic or atmospheric free surface. The intermediate coordinate surfaces 
are allowed to move as freely as one desires. Thus, external gravity waves can be rep- 
resented by the top coordinate surface, and meteorological fronts and inversion layers in 
the atmosphere and refractive layers in the seas can be represented with enhanced reso- 
lution by the internal quasi -horizontal coordinate surfaces. The action of subgrid mixing 
has been represented as true tensor quantities, which corrects a deficiency in the intuitive 
subgrid mixing formulation in that rigid-body rotation of the atmosphere under adiabatic 
conditions is an allowed solution. 

Langley Research Center 

National Aeronautics and Space Administration 
Hampton, Va. 23665 
April 1, 1975 


19 



APPENDIX A 


GENERAL TRANSPORT EQUATIONS IN COORDINATE SYSTEMS 
WITH TIME -DEPENDENT METRIC TENSORS 


In this appendix, the conservation equations for continuous media are derived in 
tensor notation. The control volume is allowed to move relative to the rotating Earth. 

The resulting general transport equations are applicable to any conservable quantity of a 
circulation system. 

Since the control volume may be in motion with respect to the Earth, the Earth rela- 
tive velocities u^ are related to the control -volume relative velocities v^ by 

u i = v i + s i (Al) 

where s^ denotes the velocity of the control volume relative to the rotating Earth. 

Consider a conservable intensive quantity Q inside a small coordinate volume v. 
The total storage rate of Q inside v is given by 


_d_ 

at 



(A2) 


where p is the density of the fluid in v and 
The outward flux of Q across the surface S 


a is the determinant of the metric tensor, 
of v is given by the surface integral 



n j dS 


(A3) 


where nj is the unit vector normal to dS positive outward. By the generalized Gauss' 
theorem of reference 10, the flux of (pQ) across S is equal to the integrated divergence 
of (pQvj) in v ; that is, 



n. dS = 


f (pQvi)|j^-d„ 


(A4) 


Conservation of Q requires that the storage rate plus the transport rate equal the rate 
of creation of Q in v, or 


(pQ) >/ adi' + J (pQv^)| j^a dv = J S'(pQ/a) &v 


(A 5) 
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where S’ is the creation rate of pQ^a" in dr. If the requirement that the control 
volume be constant in the coordinate space (but not in the physical space) is imposed, then 
the time differentiation of the first term in equation (A5) can be carried inside the volume 
integral. Since the resultant equation must hold for an arbitrary v, the integration can 
be dropped to give 


^ (pQ^a) + (pQv^j j/a = S’(pQ^a) 


(A6) 


Equation (A6) may be rewritten as 


^ (pQ^a) + - ~r (pQv j y/a) = S’(pQ/a) 
9a 


(A7) 


Equation (A7) is the general transport equation for a conservable intensive scalar quan- 
tity Q. 

The derivation of the conservation equation for tensors yields additional 

terms. In the present report, such terms are grouped with the source term to form a 
pseudosource term S^apZ^' ’ ‘ j . Thus, a common form for all conservation equations 
is obtained, which is 


it (^ 4 : : :) +^(^ 4 : : : ^ = : :) 


(A8) 


The derivation of the conservation equation for momentum is presented in the fol- 
lowing discussion. Newton's second law for the infinitesimal control volume dr, if con- 
sidering only the force of the pressure gradient and gravity and neglecting viscous forces, 
is 


du. 

p ^ d ldT + ( a c)J 


= Va dr(-P ( . +g.p) 


(A9) 


where (a c ) . is the coriolis and centripetal acceleration arising from the Earth's rota- 
tion and duj/dt is the comoving time derivation of uj, which, from reference 11, may 
be expanded in terms of the velocity relative to the control volume s 1 = u 1 - v 1 as 


dUj 

dt 


du 

IT * v U) 


i|J - s |i u j 


(A10) 
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Thus, equation (A9) can be rewritten as 


di/\/ap 



+ v J u. 


i|j 



= Va dt'|-P| 


E- 


u 


+ Sjp 



(All) 


Equation (All) can be rewritten in conservative form by the following procedure. 
First, in equation (A6), Q is set to unity. The resulting equation is the mass- 
conservation equation 



(A 12) 


Next, equation (A12) is multiplied by u^ and added to equation (All) to yield the conser- 
vative form: 



Equation (A13) can then be rewritten in the form of equation (A 8) by the following 

procedure. First, the term ^pu^v^ is expanded by identities to yield 

1 3 

^V 3 ), . = ^(v^PV 3 ) ‘ (i m j) (v/^uj (A14) 

If the following terms are collected in the pseudosource term 


S (v^p«i) = \/a£p u + gi p - P ( a c ) i + P s ]i u jj + (\/apv j u m ) (A15) 

then the momentum -conservation equation (A13) takes the desired form 



(A16) 


which is an alternate form of momentum conservation expressed in equation (A13). 
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Since equation (A16) holds for the infinitesimal volume Av, it can be integrated 
over a finite volume v. The result must still hold, although the resulting equation does 
not represent a true tensor. In numerical analysis, the solution to a given problem is 
often obtained by breaking the problem space into finite volumes. Equation (A5) for 
scalar transport, for example, might be solved by such a technique. The numerical 
approximation of equation (A5) causes some information (fine structure) about the state 
of the system to be lost or ignored. Unfortunately, the fine structure can have a signifi- 
cant role in the solution. These intuitive ideas can be symbolized by the relationship 

Q = <Q> + Q' (A17) 

where ( Q) is the average value of Q taken over the gridcell volume and time step, 
and, thus, Q' is the difference between an instantaneous point value of Q and the aver- 
age of Q over the four -dimensional space-time gridcell "volume." A model employing 
a discrete space-time grid represents the field of Q by estimating (Q). Although the 
space -time gridcell average of subgrid variations Q' vanishes (as can be seen by 
applying the operator ( ^ to eq. (A17)), subgrid variations must still be accounted for 
in estimating as shown in the following example. 

Applying the operator ( ) to equation (A6) and substituting the generic expres- 
sion (A17) for Q, vJ, Q|j, and vj yields 

(Jt^^)) + (< Q > +Q ') + v], y) 

+ <^<v/Sp)[<Q|j) + (Q|j)’](<'' ) > + v J f) 

+ <^</Sp)(<Q> + ('{))]) (/Spq)) (A 18) 

where subgrid variations of p and & are neglected. By definition of the 
operator ( ) , (Q) is constant for general Q throughout the space-time gridcell; 

thus, equation (A18) becomes 
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(^(pQV5')) + ^/ip lj )[<Q><vi> + <Q><vj’> + <vi><Q'> + <Q'vl'>] 
♦ </5p>[<Q |j )<vO +< Q| j) <vJ , > + <vO((Q |j y) t <(Q |j )V)] 

+ (-Sap) (Q) + (Q) ^( T jj) ^ <Q’> + (q’ (vjj) y 

= (s-(/Spq)) (, 


Applying the operator ( ) to equation (A17) yields 


(Q) - «Q» +(Q'.) - <Q> +{Q') 


Thus, 

<Q’> =0 (A20) 

with the result that equation (A19) reduces to 

(£( pQ ^)> + <v^ p ii)[<Q> < vi > ♦ <Q'v r >] ♦ <^)[<Q|j) <v J > 

+ (( Q li) vI '3 <Q> ( V !j) + ( Q ’( V lj)) -<S'(^pQ)) (A21) 

The operator < ) is an integral operator; therefore, for general Q, 

= (A22a) 

(by neglecting variations in the Christoffel symbols of the second kind over a space-time 

gridcell). Similarly, by use of equation (A17), 

( Q ij)’ =Q ir< Q u> < A22b > 
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and 


(QV(«-< Q >)u=V< Q >ii 


or, by use of equation (A22a), 

(Q, >ir Q ir< Q b> 

Thus, by use of equation (A22b), 



Applying equations (A22a) and (A22c) to equation (A21) yields 


(^(pQ/a)) + /a(<P> <Q> <v j >) 

I J 

= -^((p) (q’v 1 ’)) . +<^S’(/apQ)^) 

f J 


(A22c) 


(A23) 


The first right-hand term of equation (A23) can be interpreted as the contribution to the 
flux of Q into the unit coordinate volume of the statistical correlation of the space-time 
fluctuations of Q and v^. The modeling of this term is described in appendix D. 

Setting Q = Uj in equation (A23) and expressing the source term S' by equa- 
tion (A13) yields the space-time averaged form of the momentum -conservation equation; 
thus, 


* /a(< p) < «i> < V J >) ( . 

- -/s(<P> <"l'vJ’)) +^S , (/apu i )) 

1 J 

= -Va (<P> + «n P - P(a c ) i + P s |i u j] ) 

where s^ = u^ - v^ is neglected and v^ is replaced by u^ . 


(A24) 
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For convenience, in equations applied to a discretized space-time grid, the aver- 
aging brackets ( ) enclosing unprimed quantities are deleted, which should cause no 
misinterpretation because only averaged quantities are available to the finite difference 
treatment. 
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APPENDIX B 

AN APPROXIMATE METRIC TENSOR FOR NEARLY 
SPHERICAL POLAR COORDINATE SYSTEMS 

The coordinate system adopted is a variant of the standard co-latitude-longitude- 
radial system. At the land-to-air interface and the sea-floor-to-water interface, x* is 
the co-latitude angle, x^ is the longitude angle, and x^ is defined as zero. That is, 
the x^ = 0 surface follows the land and sea-bottom terrain. The x^ coordinate lines 
(which are the loci of constant x* and x^ are locally normal to the x^ = 0 surface, 
and, in general, the x^ coordinate lines are locally normal to a surface of constant x^. 
The loci of constant co-latitude and longitude, on the other hand, are radial, and so x^ 
and x^ generally deviate somewhat from the local co-latitude and longitude for non- 
zero x^. The following sketch is a simplified (two-dimensional) illustration of the 
coordinate system: 


Co -latitude, longitude constant 
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Let the local radial deviation of the x^ = 0 surface from a true spherical surface 
be denoted by 6(x-*,x^) . Then, for a displacement dx®, the first-order change in ele- 
vation (measured from a true sphere) of the x^ = 0 surface is given by 


d6 = dx 1 + dx 2 
ax 1 ax^ 

The differential arc length squared for dx^ = 0 is given by 



The x^ coordinate line is normal to the local Earth's surface; thus, for general dx^, 



The symmetric covariant metric tensor, defined by 
(ds)^ = g.j dx* dx-* 
becomes, at the Earth's surface, 



0 g 33 
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Similarly, the metric tensor can be written for positions not necessarily on the Earth's 
surface in terms of spatial derivatives of the altitude h above the ideal spherical Earth, 
where h reduces to 6 at the Earth's surface. Thus, 






(B3) 


The formulation of the metric tensor by use of equation (B3) allows the "horizontal" 
coordinate lines to be positioned or moved vertically without restriction. 

The square of the line element (which completely determines the differential 
geometry) becomes 


(ds) 2 = g^j dx* dx-i 




dx^ dx* + 


l r sin x 



dx 2 dx 2 


+ ggg dx 2 dx 2 + 2 


3h ah 


ax* ax* 


dx 1 dx 2 


Collecting terms in ah/ax a gives 


(ds) 2 = (r) 2 dx'*’ dx* - + (r sin x* - ) dx 2 dx 2 + ggg dx 2 dx 

2 


+ |^ r dx 1 + ^dx 2 
, ax 1 ax z 


(B4) 


where the first three right-hand terms are the square of the spherical polar line element, 
and the last right-hand term is the square of the height increment along the surface of 
constant x^. Thus, if the square of the horizontal increment (the first two right-hand 
terms in eq. (B4)) is much larger than the square of the height increment along the surface 
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of constant x 2 , then the square of the line element is approximated by the square of the 
spherical polar line element, and the metric tensor reduces to the spherical polar metric 
tensor 


g ij S 


(r)' 

0 

0 


t r sin x*) 
0 


0 

0 

g 33 


(B5) 


for 


ih dxl + ih dx 2 l 


, 3x 


dx* 


« (r) 2 dx* 


dx 


♦ ( 


r sm x 


dx 2 dx 2 


(B6) 


The condition (B6) for the approximation of the metric tensor by equation (B5) is not very 
restrictive because, in a discrete grid model, the slope terms 3h/3x a are character- 
istic of a coordinate cell, rather than of a single point. For example, if 6 is the maxi- 
mum change in mean elevation at the x 2 = 0 surface from one coordinate cell to an 
adjacent cell, then the condition (B6) for equation (B5) to hold (at the x^ = 0 surface) 
becomes 


(s) 

(5) 


2 

2 





J 


(B7) 


1 2 

where Ax and Ax are the cell horizontal coordinate dimensions in radians. Thus, 
the smaller the grid physical increments, r Ax* and r sin x* Ax 2 , the more restrictive 
conditions (B7) become. More generally, if h is the maximum change in mean elevation 
of a constant x 2 surface from one coordinate cell to an adjacent cell, then the condi- 
tion (B6) for equation (B5) to hold becomes 



(B8) 


Generally, it is desirable and computationally efficient to maintain a uniform horizontal 
spatial resolution of all variables. The usual method of achieving uniform horizontal 
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spatial resolution is to filter spatially all variables in the horizontal at half -wavelengths 
shorter than the maximum horizontal grid physical increments. An effect of the spatial 
(low -pass) filter is to make conditions (B8) much easier to satisfy than would be the case 
without filtering, especially near the poles, where sin is small. Consider a complex 
Fourier expansion of h around the Earth along a parallel: 


h = £ h n= ^ An(t)e n 


i^r sin x 1 x 2 
A. 


(B9) 


n=0 


n=0 


where X n is the physical wavelength associated with the nth Fourier term and A n (t) is 
a complex function of time. For simplicity, a perfect low -pass filter is assumed which 
truncates the series of equation (B9) such that h n = 0 for A n < 2r Ax 2 (twice the grid 
physical increment at the equator). The filtered h, h', becomes 


m m -i — 

h ’ = £ h n = ^ A n(t) e n 


-i^ r sin x* x 2 

A v 


(BIO) 


n=0 


n=0 


2 

where X m is the smallest X n greater than or equal to 2r Ax' . The maximum change 
in |h’| over a grid increment |Ah’| is realized if J A m (t)| = i 5h and all other 
A n (t) = 0, where 6h = (Maximum h) - (Minimum h). Thus, 




/ -illr sinxlx 2 \ 

Ah' 

< 

\ *She Xm J 

Ax^ 


ax 2 \2 / 


1 -i r sin xn 

< 

1 fib r sin vl 

1 “J. — J. Oi.ll A I 

\ x m / 


UIl X &X11 A 

^ 2r Ax z 


or 


| Ah' | s - 1 5h 7T sin x*| 


(Bll) 


Replacing h in conditions (B8) by |Ah' | yields the new condition 
(Ah')^ « (r sin x^j (ax^) 


or, by use of equation (Bll), 
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^(6h) 2 « (r) 2 (Ax 2 ) 2 


(B12) 


The first inequality in conditions (B8) is retained: 



« 



2 


(B13) 


Therefore, the inequalities (B12) and (B13) are sufficient conditions for equation (B5) for 
the metric tensor to hold when Fourier series terms for h of physical half-wavelenths 
shorter than the east-west equatorial physical grid increment are deleted. If Ax* = Ax 2 , 
which is the usual case, then condition (B12) implies condition (B13), because 1 6h| a| h|. 

Consider the extreme case where 6h = 8 kilometers, which is certainly sufficient 
for representing, by surfaces of constant x 2 , terrain and sea-bottom topography, meteo- 
rological fronts, and refractive layers in the seas. If the left-hand side of inequality (B12) 
is required to be 1 percent of the right-hand side, then the minimum grid increment Ax 2 
consistent with inequality (B12) is 


Ax 2 = - — (10) radians = - — — (10) radians 
2 r 2 6400 


or 


Ax 2 = 57.3 - (10°) = 1.12° 

2 6400 


which corresponds to a very fine grid. 

If condition (B12) is written in equation form, then 

e (p 4“ ^ 6h ) 2 = ( r ) 2 ( Ax2 ) 2 

where the error fraction in the square of the line element is «1; then, the maxi- 
mum | 6h| is 

1 6h| (kilometers) = ^ T ^ 80 — — (^ e S^ ees ) = 71.4 Ax 2 (degrees) (B14) 
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TRANSPORT EQUATIONS FOR GEOPHYSICAL 
CIRCULATION PROBLEMS 

In appendix B and approximate metric tensor was presented for nearly spherical 
polar coordinate systems. In this appendix the simplications resulting from that metric 
tensor are applied to the transport equations developed in appendix A. 

In principle, the metric tensor allows equations which are identical for all coordi- 
nate systems. The physical assumptions and solution techniques used by the analyst 
reduce the generality of the governing transport equations to special cases, and thereby 
introduce real differences into the transport equations for geophysical circulation prob- 
lems. The specialization to atmospheric- and oceanic -circulation problems will be 
made so that one set of transport equations can describe both systems. 

The prime difference, for modeling purposes, between the oceans and atmosphere 
lies in the different density characteristics. The atmosphere is compressible, whereas 
the oceans are essentially incompressible. However, if external gravity waves in the 
oceans are considered, then from the viewpoint of mass in a given fraction of a vertical 
column, the oceans and atmosphere are similar. Consider, then, a new variable £ to 
replace p so that the similarity of the fluid bodies is emphasized. A good candidate 

for £ is the grouping P\Jg%% because p and usually appear together in the 

transport equations, and £ (generalized density) has similar behavior in the atmosphere 
and oceans. 

From equation (B3), the determinate of the metric tensor is 



where A is the determinate of the metric tensor for the quasi-horizontal coordinates. 

Applying the definition £ = p^g 33 and equation (Cl) to equations (A12), (A23), and (A24) 
yields 

^ C ) + (/fe Cvi ) li ^ =0 (C2> 


33 



APPENDIX C 



d_ 

at 



1 

\/^33 



^ g 33 ^ 



g. JL_ 


. _L_ 

\/ g 33 




S 



(C3) 


(C4) 


Equations (C2) and (C3) may be rewritten, respectively, as 


< C5 > 


and 


51 «) ♦ **>) - s ’(^ «) ' < C6 > 


■j T -j f 

where v J was replaced by u J as in equation (A24). Expanding 


and 


in equation (C4) by identities similar to the procedure leading to equation (A14) yields 
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it ?u i) - ?Uivl ) = § ?o 0 ■ ? < u i’ ui ')) <C7) 


where 


S (y/A Cu,) = y/A P (i + <g, - C(a c ) t + CsjjUj + ^, m jj (kXJ) 

—I 

+ {r} (K" f v>) 


(C8) 


By inspection of equations (C6), (C7), and (C8), it is evident that equations (C5) 
and (C7) can be obtained by replacing Q in equation (C6) with 1 and u i} respectively, 

and by computing S^\/a £u^ by equation (C8). Thus, it is allowable to speak of one 
transport equation (eq. (C6)) for all conservable quantities in the atmosphere or ocean. 

From equation (B5), an approximation to the metric tensor is 


g ij ~ 


(r r o 

i 

0 (r sin x*) 


0 

0 

^33 


(C9) 


Thus, the determinant of the metric tensor is, approximately. 


— Q r ) 2 sin x3 g 33 


(CIO) 


and the determinate of the metric tensor for the quasi-horizontal coordinates is, 
approximately, 


A = Qr ) 2 sin x-|] 


(Cll) 
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The full set of Christoffel symbols of the second kind are, by symmetry in the lower 
indices, 


{/3 a y} {Vr} G“s} G 3 a} G^s} G 3 3^ 

o 

In reference 10, page 67, it is shown that, for dx locally normal to the surface of 
constant x 3 , 


w 

1 8g /3y 

2g 33 ax 3 

W 

- 1 e ay 9gy ^ 
2 ax 3 

W 

1 ag 33 
2g 33 dx a 

{3 “3} 

1 ay 9g 33 
"“2 s 3x y 

{a 3 3} 

1 9g 33 

2g 33 9x 3 


(C12) 


For the metric tensor given by equation (C9), dg pyj dx ^ = 0 and S ay = 0 for « 4 Y, 
and g aa = l/g aa \ thus, equation (C12) becomes 
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(A) 

{.“ 3 } 

(A) 

w 

(A) 


= 0 
= 0 

1 9g 33 

" 2 ^33 9xQ! 

1 8 ^33 

2g aQ! dx a 

1 8g 33 

2 g 3 3 8x 3 






(C13) 


Thus, for the metric tensor given by equation (C9), the nonvanishing Christoff el symbols 
are contained in the set 


(/3 “y} {? 3 <*} {3 “3} (3 3 3} 

The first subset in set (C14), f ^ 01 , expands to 



by use of the relationships 

g“* = 0 (j 4 a) 

aa _ 1 

® g 

& Q!Qf 


(C14) 


(C15) 


(C16) 
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Of the derivatives of the metric elements in equation (Cl 5), only 9g ^2 does not 
vanish (in general) by use of equation (C9); thus, from equations (C9) and (Cl 5), the non- 

,ing (A) 


vanishing 


are 


sm x 1 cos x 


t* 2 ) 


Then, the nonvanishing Christoff el symbols from equations (Cl 3) and (C17) are 

^2*2} = " sin x * cos x * 

(1%} - (2 ^ 1} = cot 

f 3 A _ r 3^ 1 ag 33 

l 3 a J l a V 2g„ 


(C17) 


dx a 


(CIS) 


1 9g 33 


2g d X a 


{.*.}■■ 


2g 33 9x 3 




From equation (B4) for the square of the line element (ds) , 


^-5 


(C19) 


The spherical polar system is obtained (from eq. (B4)) by setting dh/dx a = 0, for which 
case the surfaces of constant x 3 are concentric spheres. The vertical distance between 
any two such surfaces is, therefore, constant; hence, 9s/9x 3 is constant over a surface 
of constant x 3 , and, thus, 


= 


9x' 


a 


d / ds\ 
dx a \9 X 3 j 


= 0 
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or 


!f3=Uo 

dX 01 


(C20) 


for a spherical polar system. By use of equations (C18), the nonvanishing Christoffel 
symbols reduce to 



(C21) 


for a spherical polar system. 

The Christoffel symbols are used herein only in expressing covariant derivatives 
of contravariant vectors and covariant vectors; the conditions under which the spheri- 
cal polar Christoffel symbols of equations (C21) can replace the nonvanishing set of 
Christoffel symbols of equations (C18) are examined in the following discussion by com- 
paring terms in the covariant derivatives of vectors. The nonvanishing set of Christoffel 

symbols of equations (C18) reduces to the spherical polar set if and 

or equivalently 8 gg 3 /' 8 x Q! , can be neglected. From equations (C18), 



3 *\ 1 9g 33 

V " 2 g 33 ax « 


(C22) 


(3 


a 


3} 


1 8g 33 

2 g Qx a 
6 Q!Qf ° X 


From equation (C19), 



(C23) 
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or 


9 / i— — \ _ 1 8g 33 _ 9 / 9s \ 

W 33 / ' 2^ dx a 9x 3 \9x Q! / 


Thus, 


9 / 9s 


9g 33 _ g /£ — 

9x a ’ 33 9x 3 lax® 


(C24) 


By use of equations (C22) and (C23), 


f 3 "l _1 9 / 9S 

l 3 “/ "^9x3^ 


(C25) 


cn \/ g 33 9 


(•■*) 


9S 

g aQ! 9x 3 l9x Q! 


(C26) 


The line element ds, from equation (B4), is 


ds = 


(r) 2 dx* dx* + (r sin x*) dx 3 dx 2 + ggg dx 3 dx 3 

1/2 


+ (J* dxl + * dx 2> 
\dx l 9X Z 


(C27) 


Thus, 


9s 

1 


dx 


(r ) 2 ♦ 

9x J ' 


1/2 


(C28) 


9S 

9x 2 


(rsinxl) +{5 


1/2 


(C29) 
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and 



or, by use of equation (B3) for the metric tensor, 



Thus, equations (C25) and (C26) become, respectively, 




Consider the covariant derivatives of a general contravariant vector 
corresponding covariant vector Tj where 


T* — 

T |j 


8T 1 

9x3 


+ T 



which is equivalent to 


a _ 9T 
1 3 a x 3 





(C30) 

(C31) 

(C32) 

(C33) 

(C34) 
1 and of the 

(C35) 
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9x^ 



Also, 


T 


9T; 


‘li axi 


which is equivalent to 





T 


sir 


9Tg 

axJ 



From equations (C18), (C35), and (C37), 


(C36) 


(C37) 


(C38) 


a 8T a 
I £ 9x^ 



(C39) 


aT 


a 


a 

I0%x0 



(C40) 


Thus, and do not depend on 

spherical polar coordinates. 

From equations (C35) and (C18), 




and are the same as in 




8T a 

ax 3 


+ T" 



(C41) 
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or, from equation (C34), 


n/^33 


,a _ 8 T a _ g aa 8xC 
1 3 9x 3 


9h 


ax 


a 


'\Js 


act 


Multiplying equation (C42) by yjs aa ^js^ Y ields 


(C42) 


\/ g 33 


13 


\/ s aa 9T“ 1 

t 3 \[% 33 

\/g 33 9x3 2g aa ^33 

df a is 8 

1 ah \ 2 


~ “ o 1 3 ~ 

8Zg Z 9Zg 

k)J 



ax c 


ah 


ax' 


a 


(C43) 


where and dz^ are the physical components of T-* ^or and dx-*, respectively, 

and the approximate metric tensor given by equation (C9) is assumed so that ^i = ^ii T 

and ”3 ~ The second right-hand term in equation (C43) is the contribution of 

thus, in equation (C41) for t“ 3 , can be neglected if the square of the 

physical slope of the constant x 3 surface (a h/az^j 2 is negligible; or, from equa- 
tion (C41) ^to the order (dh/dz^ ^ 

T ?3 = S +X3 {^?} 


8 T 

dx c 


a 


(C44) 


which is also the order of approximation of the metric tensor by use of equation (C9), as 
discussed in appendix B. In the remainder of this appendix, the symbol = denotes 
approximations accurate to the order f ah^/az^ 3 . 
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From equations (C36) and (C18), 



or 



where, again, 



(C45) 


(C46) 


(C47) 



(C48) 


with equation (C48) being exact for i = 3. Expanding the second right-hand term in 
equation (C46) gives 


T » / * \ If 1 8g 33 

3 [fo) ' 2 3 g3 3 3/2 «■“ 



by use of equation (C22) for 



T 3 ~ 1 9 ^ 3 


Therefore, equation (C46) becomes 


(C49) 


which is equation (C45) for negligible 
and (C33), 



and 9g 


33 / 9xQ! - 


From equations (C36) 
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_3 _ 9T ,py 

T IS 9*3 +T 


{Vs} +t3 (3 3 3} 


aT° 

9x 3 


+ T 7 


\J % 33 


3x J 



+ T 0 



(C50) 


or, by using equations (C47), (C48), and (C18), 



Thus, by use of equation (C50), 



to the order of ^ah/az^^; and, furthermore, by use of equations (C51) and (C47), 


(C51) 


(C52) 
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T 3 ~ aT 3 - 1 8T3 

^ 3 dZn v/goo 9X 3 




(C53) 


to the order of [dihjdz^ . 

In summary, in the computation of the covariant derivatives of contravariant vec- 
tors T}., aggg/ax®, 1 3 3 \ , and { 3 ^ 3} can be neglected with errors no greater 

* 3 ~ \2 
than the order of the square of the physical slope of the constant x° surface (9h/3 z 1 . 

By use of equation (C40), is not dependent on ^3 3 y| > or { 3 y 3 } • 

Now, in considering Tgj^, T^g, and Tgjg, 


-1 rri^ p. ~ rpl ~ 

i|j |j s ik |j g ii 


(C54) 


by use of equation (C9) for the approximate metric tensor. Then, by use of equations (C49) 
and (C54), 


g 9Tg 

T 3|q! ~ T | a g 33 - \[%33 


8x 


(C55) 


According to equations (C38) and (Cl 8), 


* 3 i 


(C56) 


or, by equation (C48), 


3 |a dx a 


{^33 ^3) “ T 3 (3 3 ck} 


(C57) 


which is equation (C55) for negligible 8ggg / dx01 and / 3 3 \ . Therefore, equa- 

if 1 jii. ^ 


tion (C56) can be written as 


3To 


3 I“ dx a 


(C58) 
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(C59) 


(C60) 


(C61) 
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In summary, in the computation of covariant derivatives of contravariant and 
covariant vectors, { o u ~.) and can be neglected with errors no greater 


3 ’ (3 3 a} 311(1 {3 01 3} 


than the order of the square of the physical slope of the constant surface (ah /azl . 

For the nearly spherical polar coordinate system with the exact metric tensor given 
by equation (B3), the nonnegligible Christoff el symbols are 


sin x 1 cos x 1 


{2 *2} - - 

{V2} = {Vl} =cotxl 

fs\ 1 9g 33 

1 ? 3 J ” 2 §33 9x3 


(C62) 


and 


3g 


33 


9x 


a 


^ 0 


(C63) 


for computation of covariant derivatives of contravariant or covariant vectors. 
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INTERPRETATION OF THE FINE -STRUCTURE INFORMATION 
AS A SUBGRID MIXING TENSOR 


In appendix A the effects of uncertainties in the state of the fluid on mean fluid 
properties are described. In the present appendix, the effects of mean fluid properties 
on uncertainties in the state of the fluid are modeled. 

The approximation of the transport equations by finite difference equations neglects 
some fine -structure information about the fluid state. In an effort to recover the mutual 
coupling between the uncertain fine -structure fluid properties and the averaged (over a 
space-time gridcell) fluid properties, the pertinent fine-structure properties are 
expressed in terms of mean fluid properties by means of a tensor formulation of the 
standard subgrid mixing model (ref. 1). The tensor formulation of subgrid mixing cor- 
rects a deficiency in the standard treatment in allowing rigid-body rotation of the atmos- 
phere under adiabatic conditions. The subgrid correlation tensor ^u-* uj'^ in equa- 
tion (C8) will be approximated from the strain-rate tensor as in reference 3. For the 
global scale circulation, Smagorinsky has modeled the elements of in the hori- 

zontal space. The simplest tensor analogy for Smagorinsky' s model is described by 
letting the strain-rate tensor be 

V - u j|k - u k|j < D1 > 

the deviatoric strain-rate tensor be 


^jk £ jk " 3 £ P g jk ( D2 ) 

and the subgrid-velocity-correlation tensor be 

KW) = *) 2 ijijk < D3 > 

where Kq is the Von Karman constant, A is the coordinate increment, and 
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Consider next the subgrid transport term , where the tensor analogy to 

Smagorinsky's model is 


p<“i'«') - -p[#( K o A ) 2 ^|j 


The deviatoric strain-rate tensor can now be expanded explicitly. 

From equation (Dl) for the three-dimensional case, f ^ is given by 

— >i 

2u l|l ( u l|2 +u 2|l) ( u l|3 +u 3|l) 

£ jk = ( u 2 |l + u l| 2 ) 2u 2 | 2 ( u 2|3 + u 3 | 2 ) 

( u 3|l + u l| 3 ) ( u 3|2 + u 2| 3 ) 2u 3|3 

which, from equations (C62) applied to equation (D2), expands to 
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(D7f) 


(D7g) 


(D7h) 






(D7i) 


The hydrostatic assumption neglects unbalanced forces in the quasi-vertical 
direction x^. In order to specialize equations (D7) and (D3) to the quasi-two -dimensional 
case consistent with the hydrostatic approximation, the following conditions are imposed: 

(1) The quasi-vertical transport of subgrid quasi -vertical momentum p^u^u^’^, 

or equivalently is neglected. 

(2) The quasi -horizontal shear of quasi -vertical velocity is negligible compared 

with the quasi-vertical shear of quasi-horizontal velocity, or 
Equations (D7) for £ j become 
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^ + -^( U 2s/S22 ) ( sin x 1 ) 2 - 2 sin x 1 cos x 1 |-^L] 
V g 11 g 22 ' 'j 1 W g 22J 


sin x- 1 


\[%22 


Uj g + (U 2 sin x 1 ] - 2 cos x* U 2 

\ / ,1 


(D8b) 




( sin x*) 


(D8c) 


-2 *1 
£ 2 = - £ .l 


(D8d) 


■ 3 ~Ki 


(D8e) 


i® 


.1 


’33/ 


.3 


(D8f) 


i 2 s = ^| 

/S 22 


(D8g) 


£3 _/!22U2 

.2 goo -3 

\ & 33 / 


(D8h) 


The quantities D T and D g from reference 3 are related to and by 


and 


+ 2 cos x 


~ 1 
r ± 


l-S-.o 


\[%22 


f 2 1 U 2 

— + 2 cos x 1 t = Dg 
sin x 1 ^22 " 


(D9) 
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In Smagorinsky* s subgrid mixing model, the second left-hand terms in equations (D9) 
cause rigid-body rotation stresses, contrary to physical principles. The second scalar 

of £ l . is taken in the tangent plane in equation (D3) for the quasi-two-dimensional case. 

* J 

Thus, 


H = 



1 

.1 









(DIO) 


For rigid-body rotation west to east, Uj = 0, U2 = c sin x-*-, and Ug = 0, where 
c is a constant. Substitution of these velocities into equations (D7) for the three- 
dimensional case and into equations (D8) for the quasi-two-dimensional case results in, 

respectively, f 1 - = 0 and l a R = 0; thus, 1 and the subgrid velocity-correlation (or 
Reynolds stress) tensor vanish, by using equations (D3), (D4), and (DIO) for both cases. 
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APPLICATION OF TENSOR TRANSPORT EQUATIONS TO 
THREE-DIMENSIONAL MOMENTUM AND SCALAR CONSERVATION 
EQUATIONS IN PHYSICAL COMPONENTS 


In this appendix the tensor transport equations and associated source terms are 
converted into the more familiar physical components. The conversion begins with the 
mass transport equation (C5) which is given as follows: 


W 1 ? ) + ^ ?vI ) - 0 


(El) 


Since A is not a function of time and v J can be converted to Vj 
cation by equation (El) simplifies to 


merely by multipli- 


9£ + JL 

sfA 9x3 




= 0 


(E2) 


The next transport equation from appendix C to be considered is the internal-energy 
equation. The first law of thermodynamics states that the increase in internal energy of 
a system results from heat being added directly to the system and/or from work being 
done on the system. If the heat added to a unit of mass m is Q, and the work done 


on m is the pres sure -volume work of compression 
tion for internal energy is, from equation (C6), 


.PJ. 

p V 


then the conservation equa- 


£(n/A ?e) + ^-j(\/a ?Ev j ) = s(\/A ?e) - JL (>/A ?<E’J’>) (E3) 

where 

S (\/A ?e) = \/A ?Q - \/Ayi^ Puj j (E4) 


and where Q includes radiation heating as well as latent heat release from condensation 
of water. 
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The conversion of equations (E3) and (E4) consists of dividing by \Ja, expanding ujj 
to -p -^-r (^a u-*), and replacing vi and by Vjy^gTT and respectively, 


to get 


where 


£ (?e) (/I ?e) - k s (f 




(E5) 


For the atmosphere, the specific humidity conservation equation has rainfall as the 
negative of its source term S (,/A from equation (C6), 


*" * " &)- , <' K “> 

Consider next, the momentum equations (C7) and (C8) 

S^PH) P-H) = S(^A C U.) - • (^K-i'uO) 


(E6) 


(E7) 


and 


S (n/a CUj) = \/A P u + ?gj - ?(a c ). + ? S ;. Uj - ?( vJu m + ( uJ ’v)) 


(E8) 


From equations (C62), the decomposition of the Christoffel symbols into the hori- 
zontal and vertical spaces yields 


&■>}■{.■} (•••) 


(E9) 
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where only g 33 in is a function of time. The term \/A CsJ.Uj inequa- 

tion (E 8 ) expands to 

+ ^"» +, ' , G 3 ‘}" 3— 3 G 3 J - 3 ) 

"^?|j u 3 * s 3 { 3 3 i}vi 3 ) (E10 > 


, (3 a ij ~ °» and {3 3 i} “ 0 f or 1 ^ 3 ( by use 

of eqs. (C62)). For the U 3 momentum component, the first left-hand term in equa- 


1 2 

because s = s = 0 everywhere 


tion (E7) becomes 


{'fc ? u 3 ) ?U 3\/®33 ) 


3 ) + ^ ? u s £(1/^) 

= Tt ?U s) + ^ ?U 3#S3 s l 3 


which follows from defining the physical increment 




(Ell) 
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s 3 
S | 3 


(El2) 


Expanding further, by use of equation (E10), 


^s) = 1^33 ^('/A CU 3 ) + l/A u 3 + 4%} u 3 

■\Z S M5f(^'' U s) + '^ !: (^5 U 3 + 83 {> S 3} U = 


(El 3) 


Then, the third momentum equation in physical components is, from equations (E7) to 
(E10) and (E13), 



(E14) 


The second term on the left-hand side of equation (E14) expands to give 


57 



APPENDIX E 


I'/ac^U:- L - 






+ /a ?u, 


8 


»JJ 2 ft™ 8x1 


j( g 33> 


(El 5) 


From equation (C63), the first two terms of the second summation on the right-hand side 
of equation (El 5) are negligible; thus, from equations (C62), 


^ ?U 3 A^ (g33)= ^ U3Vs ^ 


(E16) 


Thus, the right-hand term in equation (El 4) 


^(3 3 3 } ?( 


V 3 U s) 


can be canceled with the second term of the expansion (El 5) of the second left-hand 
term of equation (E14). Similarly, cancellations can be made for the following term in 
equation (E14): 


\/A { S 3 3 } ?(< U 3 ,U 3’)' 


and a portion of an expansion (following the format of eq. (El 5)) of the right-hand term 


to leave, after division of equation (E14) by \Jg 33 /a, 


ifcu,) + — (/a ru, = -Pi, + C - c __5 - -.-(/a t:/ ij ' i- 

’ ^ 8x, \ 3 ^ij 13 m ^ ^ ax, \ \ >1* 


(El 7) 
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Thus, 


U(/a C u 3 ) = 


Va 


S3 

-P| 3 + C^-5 


( a c) 


3 19 


^33 \P 33 ^ 9xj 


T Vac/ u 3 


u/ 


' 3 


(E18) 


For the momentum component along dx 1 , ^|g 1;l = R is constant. Substituting 

Ui = into equations (E7) and (E8) yields, by use of equations (E10) and (C63), 




(E19) 


and 


( a c) 1 


T"3 

Ja \ > v/elT ® x fin fin fin ax 


io 9So 

+ — s. 05 




of oi a 




_1 9_ 

\[A 9xj 


/3/3 


1 /3 


*( V 0 V « + <V D c'>) 


/__ / u,’ 

Va ?/ Uj' -J 

\ >/* 


(E20) 


Equation (E20) can be written more explicitly by using the result of appendix C; 
that is, of A 01 &) ’ ° nly 2 0 = COt Xl 1S nonzer °‘ Thus > ec l uation (E20) reduces to 


fA 


/ \ \j goo 3P ?§i C ) 1 P 9 ^3 

! [Ja SU,) = -3L-3 3 JP + i i + _S f U„ 

' ' >/^i i ax ^ ^ 


- 1 E </ J=)) + ^ cot X ‘( V2U2 + < U2 ' U2 ’» <E21> 
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Conversion of the momentum component along dx^ is accomplished in a similar 
manner, except that now extra terms are obtained from the partial derivatives in the 

horizontal direction because Ug = u 2^/&22 and \]&22 ~ r sin x ** When the substitution 
for U 2 has been made in equations (E7) and (E8) and the results converted again to the 
form of equation (E7), there follows, from equations (E10) and (C63), 



(E22) 


and 


J-sfi/AtU ) - 8F 1 ?g2 

. 'i^aa f a 


/ g $3^ g 22 ^ 


{*%} ?(V« + <W» 


u; 




\p22\lH1 


5( U 2 V 1 + < U 2' U l'» 


(E23) 


Equation (E23) can be simplified further by expanding 01 ^ 
appendix C that 


and noticing from 


^2 ^ 2 } = “ sin x * cos x * 


and 


{ 2 2 i} = cotxl 
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are the nonvanishing terms. When these Christoffel symbols are substituted into equa- 


tion (E23), the term containing 
is 



cancels the last term and the simplified result 


-L sux = . £SE J* + if?. . + _jl_ !?§ u 




-^cotx 1 (v 1 U 2+ (u 1 'U 2 ')) 


(E24) 


This completes the set of conservation equations governing the motion of a fluid 
with three momentum components. The system allows for the time dependence of 

^ggg, which can be computed from 33^ = \j 633 ^3 by equation (Ell) if s^ is 

specified. 
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RELATIONSHIP BETWEEN THREE-DIMENSIONAL MOMENTUM 
EQUATIONS AND HYDROSTATIC MOMENTUM EQUATIONS 

When the full set of transport equations derived in appendix E are integrated on a 
digital computer to simulate the seas or atmosphere, the horizontal scale is generally 
much greater than the vertical scale. Consequently, weather prediction or ocean simula- 
tion would be prohibitively expensive if approximations were not made to alleviate the 
short time steps associated with integrating the vertical momentum equation. 

The classical method to obtain a long time step is to neglect accelerations along 
the gravity vector, computing pressure by the hydrostatic approximation and then com- 
puting vertical velocity from the energy equation while conserving internal energy. (See 
ref. 1.) The loss of similarity between the resulting numerical model and the real world 
is associated with buoyant or convective overturning. (See ref. 3.) Consequently (for 
example), a hydrostatic circulation model is expected to be a poor vehicle for simulating 
thunderstorms. 

The momentum equation derived in appendix A will be specialized to the hydro- 
static approximation in this appendix. The hydrostatic momentum equation can be 
derived from equation (A13) which is revised as follows: 


d_ 

at 


(i/apuj + - /^JpsjjUj - 




Equation (A13) has three components which span the physical three-space. The hydro- 
static approximation employed herein is that the right-hand (force) term of equation (A13) 
for the nearly vertical x^-direction is zero. Then, equation (A13) becomes 


±(/Ipu a ) + Va(pu a uJ) ( . 




pS ta u j 




Tt{^ pu s) + ^( pU 3 vj )|j 


/a 


P s' 


|3 U j 



0 


-\ 




J 


(FI) 


The hydrostatic assumption 
P |3 = p S3 


(F2) 
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and the condition that the pressure P vanish at the top (x** = X^J of the atmosphere or 
seas yield 



where the gravity -acceleration components g^ in equations (FI) and (F3) are 

Si = -g (F4) 

1 8x l 

The third component of equation (FI) is no longer needed because the x^ velocities can 
now be computed from the internal-energy conservation equation (E3). 
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